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Abstract. We present simulations of colliding supersonic stellar winds with RAM- 
SES. The collision results in a double shock structure that is subject to different instabil- 
ities. Properly modelling these instabilities requires a high enough resolution. At large 
scale, orbital motion is expected to turn the shocked zone into a spiral but we find that 
in some configurations the Kelvin-Helmholtz instability (KHI) may disrupt the spiral. 
A colliding wind structure is also expected in gamma-ray binaries composed of a mas- 
sive star and a young pulsar that emits a highly relativistic wind. Numerical simulations 
are necessary to understand the geometry of such systems and should take into account 
the relativistic nature of the pulsar wind. We implemented a second order Godunov 
method to solve the equations of relativistic hydrodynamics (RHD) with RAMSES, in- 
cluding the possibility for Adaptive Mesh Refinement (AMR). After a brief overview of 
our numerical implementation, we will present preliminary simulations of gamma-ray 
binaries. 



1. Introduction 

Massive stars posses highly supersonic winds. In binary systems, the interaction of 
two winds creates a double shock structure which geometry depends on the momentum 
flux ratio of the winds. The winds are separated by a contact discontinuity. Numerical 
simulations show that several instabilities may arise in the colliding wind regio n (|Pit- | 
tard 2009 ). For isothermal winds, the non-linear thin-shell instability ( |Vishniac||1994 f 



dominates the collision region while the KHI can strongly perturb the shocked region 
between adiabatic winds ( |Lamberts et al.|20lT| ). 

Up to now, high-resolution simulations have focused on the region close to the 
binary. At larger scale orbital motion turns the shocked structure into a spiral. Matter 
is expected to be ballistic, with a distinction between both shocked spiral arms (Parkin 



et al.||20lT] ). The exact structure of the spiral is still to be studied. A spiral structure 



is observed in WR 104, a binary composed of a Wolf-Rayet (WR) and an early-type 
star. Owing to dust emission, it displays a spiral in infrared up to a few hundred times 
the binary separation ( [Tuthill et al.||2008| ). The formation of dust in such systems is 



still poorly understood, and is likely to be facilitated in the colliding wind region. We 
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perform simulations of this system to determine the structure of the colliding wind 
region and put constrains on dust formation. 

y-ray binaries are binaries composed of a massive star and a young pulsar (see 
e.g. |Abdo & Fermi Co llaboration (2009)) with a tenuous, highly relativistic wind. The 
collision between the pulsar wind and the wind from the companion star, which re- 
sults in y-ray emission, is expected to present similarities with the collision between 
stellar winds ( |Dubus||2006| ). Still, the relativistic nature of the pulsar wind is likely to 
modify the dynamics and emission properties of large scale colliding wind region that 
can be observed in radio (see e.g. Moldon et al.| ( |2011| )). To investigate the impact 
of the relativistic nature of the pulsar wind, we developped a relativistic extension to 
the hydrodynamical code RAMSES ( |Teyssier|2 002). We explain our numerical imple- 
mentation, focusing on the AMR and present some preliminary simulations of y-ray 
binaries. 



2. Large scale structure of colliding wind binaries 

We use the hydrodynamical code RAMSES for our simulations. It is a second order 
Godunov method, that solves the equations of hydrodynamics. We focus on the adia- 
batic limit. We use AMR to increase resolution at the shocks and properly model the 
instabilities while simulating the large scale structure at reasonable computational cost. 
The winds are generated following the method by |Lemaster et al.| ( [2007| ). We include a 
passive scalar to distinguish the winds and determine mixing. 
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Figure 1 . Density (upper row) and mixing (bottom row) for 2D simulations with 
increasing velocity ratios and decreasing density ratios (from left to right). On the 
left panel both winds are identical. The length scale is the binary separation. The 
maximum resolution is set by l max = 15. 
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FigJT] shows density and mixing for winds with equal moment fluxes (M\v\ = 
M2V2) but increasing velocity ratios (and decreasing mass loss rate ratios). When both 
winds are identical (left panel), the spiral arms have the same step and same width. 
When the winds have dierent speeds, the spiral arm propagating into the faster, lower 
density wind expands, while the arm propagating into the denser and slower wind gets 
compressed. Mixing is more important in the wider arm (Lamb erts et al.|20"T2 ). When 
increasing the velocity dierence even more (vi = V2 = 20), the spiral structure is de- 
stroyed due to the KHI. Performing additional simulations, we show that an important 
density dierence can have a stabilising eect, even when winds have dierent speeds. We 
find that the step of the spiral is primarly set by v x P where P is the orbital period of 
the system and v the speed of the dominating wind, with significant corrections fromt 
he weaker wind in some cases. 

We perform a detailed study of WR 104. Our goal is to determine whether a hy- 
drodynamic model with adiabatic winds can reproduce the observed spiral structure in 
WR 104. The WR wind is very hostile to dust formation because of its low density, high 
temperature and absence of hydrogen and the colliding wind region is likely to facilitate 
it ( |Walder & Folini||2003| >. In WR 104 the wind from the WR strongly dominates the 
wind from the early-type star and the shocks form very close to the latter. A very high 
resolution is thus necessary close to the binary and is incompatible with a large scale 
3D simulation. We thus perform a small scale 3D simulation and a complementary 
large 2D simulation. The 2D simulation confirms the Archimedean spiral and indicates 
the highest densities are reached at the edges of the spiral, suggesting they could be 
the location for dust formation (Lamberts et al. 2012). The results from the 3D sim- 
ulation are given on Figj2j The temperature map suggests temperature could become 
low enough to enable dust condensation (< 6000 K) at half a turn of the spiral, which 
is compatible with observations. Mixing (middle panel) between the winds is likely 
to bring hydrogen and facilitate chemistry. Still, the density map shows that, even in 
the walls of the spiral, the density remains below the critical value for dust formation 
(~ 10 6 cm" 3 ). An adiabatic model is thus unable to account for dust formation in WR 
104. Cooling is probably present, and increases the compression ratio of the shocks. 
It might also trigger the non-linear thin shell instability, which results in small high 
density zones. 





Figure 2. Density (g cm 3 ), mixing and temperature (K) in the orbital plane of 
the 3D simulation of WR 104. The length scale is the binary separation. 
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Extension of RAMSES to relativistic hydrodynamics, application to y-ray bi- 
naries 



The equations of RHD can be written as a system of conservation equations (c = 1) : 

(1) 
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where D is the mass density, M the momentum density and E the energy density in 
the frame of the laboratory. The subscripts i, j stand for the dimensions, 8 l -j is the 
Kronecker symbol, h is the specific enthalpy, p is the proper mass density, is the fluid 
three-velocity, P is the gas pressure and y the adiabatic index. The Lorentz factor Y is 
given by 

r = -f=^ ( 2 > 

^Y^v 2 

These equations have a similar structure to the equations of hydrodynamics but are 
more complex to solve because strongly coupled to each other by the Lorentz factor 
and the enthalpy. An additional numerical constraint arises from the the fact that the 
velocity must remain subluminal. The similarity with the equations of hydrodynam- 
ics allows us to closely follow the algorithm implemented in RAMSES, performing 
localised changes. 

A first difficulty arises when passing from the conservative variables (D, At, E) T 
to the primitive variables (p, v, P) T . There is no explicit solution, and a root find- 
ing algorithm is necessary to recover the primitive variables. Care has to be taken to 
avoid numerical problems in the ultrarelativistic and non-relativistic limits (Mi gnone] 
& Mc Kinney 2007). Second-order precision is implemented in RAMSES following a 



MUSCL-Hancock method. The reconstruction of the primitive variables in space is not 
affected by special relativity but the prediction in time requires the determination of 
the Jacobian matrix of the system given by Eq. [T| In case of a superluminal velocity 
in the reconstructed state, the code switches to a first order scheme. The relativistic 
summation of velocities changes the determination of the wavespeeds and the timestep. 

In RAMSES, AMR is implemented following a tree-based method where par- 
ent cells are refined in child cells in a recursive structure ( |Kravtsov et al.|[T997] ). The 
child cells are gathered together in octs of size 2ndim. Interpolation from I to I + 1 
is done at the interface between dierent levels and when creating new refined cells. 
We perform the reconstruction on the conservative variables, and switch to first or- 
der when a non-physical state is found. A non-physical state occurs when the density 
or pressure are negative or the velocity superluminal. A physical state is guaranteed 
when (D > 0, E 2 > D 2 + M 2 ). Hydrodynamical updates are only performed at the 
highest level of refinement, and variables are computed at lower levels by averaging 
values over the child oct. Although each cell from the child oct at level I satisfies 
(D > 0, E 2 > D 2 + M 2 ), it is not necessarily true on average, over the whole oct. This 
may lead to non-physical states at level I - 1. Therefore, we choose to perform the 
restriction on the internal energy rather than the total energy. This method is currently 
implemented in RAMSES. 

Fig. [3] shows two examples of numerical tests we performed with our new rel- 
ativistic code. The left panel shows the density, transverse and parallel velocity and 
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Figure 3. Left panel : Sod test (t = 1.8).Right panel : Simulation of the propaga- 
tion of a 3D relativistic jet (T max = 7.1). From top to bottom : density at t = 20, 30, 40 
in a 3D jet starting from the left boundary of the domain 



pressure in a one-dimensional Sod test. Contrary to the classical case, transverse ve- 
locities do impact the structure of shocks in RHD. The setup is taken from (Ry u et al. 
2006| ) and constitutes a very stringent test due to the high Lorentz factor (T max = 120) 
and important transverse velocity. A very high resolution (AMR levels shown in red) is 
needed to obtain a good agreement with the analytic solution (in blue). The right panel 
shows a the results of a 3D simulation of an axisymmetric jet following the setup by 



Del Zanna & Bucciantini| ( |2002| ). Both tests show satisfactory results and indicate the 
code is ready for scientific use. 

Our aim is to model y-ray binaries. We focus on the small scale structure of 
the interaction between a stellar wind and a pulsar wind. The goal is to understand 
the impact of relativistic effects both on the structure and stability of the interaction 
region.We neglect orbital motion and focus on winds with equal moment fluxes. We 
perform preliminary simulations with various values of the momentum flux ratio and 
pulsar wind velocity. They prepare a large scale simulation to determine whether a 
stable structure is possible ( |Bosch-Ramon et aL|2012| ). Fig. 4 shows the density map 
for a simulation with the speed of the pulsar wind v p - 0.5. The KHI develops in a 
similar fashion than in the classical case. The right panel shows the positions of the 
discontinuities for simulations with increasing values for the pulsar wind speed. The 
higher the value, the more the shocks are bent towards the star. This is a relativistic 
effect due to the impact of transverse velocities on the structure of shocks. 



4. Conclusions and perspectives 

We performed high resolution simulations of colliding wind binaries at a spatial scale 
never reached before. We showed that the KHI may destroy the expected large scale 
structure. Simulations of WR 104 match well with the observed structure and indicate 
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Figure 4. Left panel : density map of a simulation with equal momentum flux, 
withvp = 0, 5, the star is on the left, the pulsar on the right. Right panel : position of 
both shocks and the contact discontinuity, in simulations with different values for the 
velocity of the pulsar wind. We have v p = 0.01 (thin dotted line), 0.1 (thick dotted 
lined), 0.5 (thick dashed line) and 0.9 (solid line). 



cooling has to be taken into account to allow dust formation in this system. To model 
y-ray binaries, we extended RAMSES to relativistic hydrodynamics. Preliminary sim- 
ulations of y-ray binaries confirm a similar structure to stellar binaries. The relativistic 
extension of RAMSES allows the use of AMR and is suited for the study of gamma-ray 
bursts, relativistic jets or pulsar wind nebulae. It will be part of the next public release. 
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